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ABSTRACT 

The fusion enhancement factor due to screening in the solar plasma is 
calculated. We use the finite temperature Green's function method and a self 
consistent mean field approximation. We reduce this to one center problems, 
because in the collision of two fusing ions, the turning point where tunneling 
may occur lies far inside the screening radius. The numerical results given by 
this method indicate that screening may be slightly weaker than that obtained 
in the most recent previous calculations. 
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1. Introduction 

There remains continued interest in the solar neutrino problem which has not been 
resolved in terms of the standard solar model and standard particle physics. Neutrinos are 
generated in thermonuclear reactions, and thermonuclear reaction rates are determined 
by solar physics (temperature and density), nuclear physics, and atomic shielding. The 
solar neutrino problem is just one of several examples where there is current interest 
in thermonuclear reaction rates; others include general stellar structure and laboratory 
plasmas. We address here the role of atomic shielding, using a finite temperature Green's 
function method. Recent calculations of screening effects in the sun have been reported by 
Gruzinov and Bahcall 1998|, who also review previous solar screening calculations. Results 



of those calculations, with which we compare our present work, exhibit a scatter of 10 



to 30% or more in the effect of screening on fusion rates. As discussed in [Gruzinov and 



Bahcall 1998|, all of these previous calculations are based either on Salpeter's weak screening 



formula, which is a classical approximation, or else make unphysical approximations in 
order to go beyond the linear regime. Gruzinov and Bahcall attempted to rectify these 
limitations by including quantum effects, treating the kinetic energy in perturbation theory. 
Out model differs from theirs in that we do not rely on perturbation theory. 

The physical model we employ is based on self-consistent Hartree plus finite 
temperature, local density approximation for exchange and correlation and assumes thermal 
equilibrium and adiabaticity in the internuclear coordinates. We solve this model making no 
further approximations or expansions, except for numerics which are carefully monitored. 
All previous calculations are essentially approximations to this model, involving e.g. 
expansions in the temperature, density or quantum corrections. [Carraro et al. 1988| also 



consider dynamic effects. Thus we believe that our approach yields results which are the 
most accurate to date. These results yield fusion enhancement factors which are typically a 
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little smaller than those recently reported ( Gruzinov and Bahcall 1998| ). 



The rate of thermonuclear fusion in plasmas is governed by barrier penetration. The 
barrier itself is dominated by the Coulomb repulsion of the fusing nuclei. Because the 
barrier potential appears in the exponent of the Gamow formula, the result is very sensitive 
to the effects of screening by electrons and positive ions in the plasma. Screening lowers 
the barrier and thus enhances the fusion rate; it is more important the greater the nuclear 
charges, and thus plays an important role in the solar neutrino spectrum. 

As shown below, the classical turning point radius which enters in the WKB integral 
for barrier penetration is very small compared with the characteristic screening lengths 
of interest. Inside this radius, the barrier potential is just the nucleus- nucleus Coulomb 
repulsion minus a constant due to screening. The constant can be interpreted simply as the 
difference in free energies of the system for either united and separated nuclei. Since both 
of these are spherically symmetric, one needs only consider one-center problems. 

We present here numerical calculations relevant to the solar core. The screening due 
to electrons in the plasma is calculated quantum mechanically by a novel Green's function 
method described in a recent paper ( |Watrous et al. 1999| ). The screening due to the ions 



is treated classically but self-consistently with the electrons. The resultant enhancement 
factors for several nuclear reactions are presented and compared with earlier results by 
other authors. 



2. Fusion rate 

We work in atomic units, e = h = rrie = fc^ = 1; the unit of atomic temperature 
is 3.159 X 10^ K. As shown by [Salpeter 1954| (here we follow [Clayton 1968| ), the fusion 
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reaction rate between species 1 and 2 is given by 

/•oo 

ri,2 = N,N2 / tljiE)viE)aiE) dE, (1) 

Jo 

where Ni and N2 are the number densities of the colhding nuchdes, ip{E) is the Maxwelhan 
probabihty that the center-of mass-energy at large separation is E; and the cross section a 
is written as a product of a penetrabihty factor P and a nuclear factor a„„c, 

a{E) = P{E)anuc{E) . (2) 

In the WKB approximation, the s-wave penetration is given by 

-2^2^ j""' [Z^Z^/r + Usc{r) - E]"''dr 
^ Jr 



- p \ 1/2 
Po{E) = ( -fj exp 



(3) 



with Eb the height of the barrier, /i the reduced mass and Usc{r) < the screening 
potential; R is the "touching" radius (the top of the barrier), and Rq is the classical turning 
radius. The integrand in Eq. (1) peaks at some characteristic energy Eq. For the solar 
interior, the corresponding classical turning radius 



Ro ^ Z,Z2/E, 



0; 



is ~ 0.005 au compared with a screening radius of ~ 0.5 au. For this reason, the screening 
potential can be taken as a constant Uq = Usc{r = 0) < within the turning radius, acting 
as an effective shift in the energy: 



ri,2 

roo 



OC 

/O 



-Uo 

Then, neglecting Uq in the limit of integration. 



N1N2 / ^{E)v{E)P{E - Uo)anuciE - Uo) dE (4) 

Jo 

roo 

/ E'/'e-''/^P{E-Uo)anuc{E-Uo)dE (5) 

Jo 

/oo 
(E' + UoY/'e-(^'+^^y^PiE')ar.UE') dE' . (6) 

-Un 



rOD 

n 2 « e-^o/^A^iA/'s / ij(E)v(E)(r(E) dE . (7) 

Jo 
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Screening thus enhances the reaction rate by the factor 

/ = e-^«/^ , it 

where Uq is given by Eq. (10) below. Since Uq < 0, f > 1, corresponding to enhancement. 



3. Free Energy and the One Center Problem 

The screening potential is just the change, brought about by the approaching ions, in 
the Helmholtz free energy 

-'^ '-^internal -^ "^j W/ 

as a function of nuclear separation. Because of the smallness of the turning radius, it is 
sufficient to consider the one-center problems and identify 

f/o = F(Zi + Z2) - F(Zi) - F{Z2) . (10) 

What we mean here by one-center problems are the calculations of F when either charges 
Zi,Z2, or both lie at the center of the plasma. In Sec. IV we calculate F{Z) explicitly 
for an independent particle model approximation, akin to the temperature dependent 
Hartree-Fock method. However, the free energy has a simple and completely general (for the 
problem at hand) interpretation, which does not depend on the mean field approximation. 
It is, by definition, the work required at a given temperature to introduce the nuclear 
charge into the plasma. The nucleus only interacts with the plasma electrostatically, so 
dF = (f){Z', r = 0)dZ', and thus the work required to assemble the charge Z at, say, r = is 

F{Z)= f%{Z',r = 0)dZ', (11) 

Jo 

where (f){Z', r) is the electrostatic potential generated in the plasma due to a nuclear charge 
Z', 

0(Z',r = O)=47ry^ r'dr'[J2 pi{r',T) - p,{r',T)\ . (12) 
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Here pe{r,T) and pj{r,T) denote the finite temperature number density (recall e = 1) of 
the electrons and various background ions charge densities, respectively. 



The one-center screening problem is solved by the methods described by |Watrous"~et 



|al. 1999| for the electrons with a Kohn-Sham formalism, using the finite temperature local 



density approximation ( [Mermin 1965| ) for exchange and correlation. Screening due to ions 
is included as well as that due to electrons. The ions are treated classically according to the 
Debye-Huckel method. Thus the electric potential generated by electrons and ions becomes 

tM^g+/#/^'^'<^';^)-^'<^'^). (13) 

r J |r — r'l 

The second term is the function denoted by (l){Z, r) above. Each ion density is given by 

p,(r-,T)=p,(oo)e-^^*W/^. (14) 

Since $(r) is positive, ions are pushed outward and electrons are drawn toward the nucleus. 
For notational simplicity in the following, the temperature dependence of the densities will 
be understood and not specified explicitly. 



4. Alternative Mean Field Derivation of the Ftee Energy 



We now demonstrate the equivalence of the F defined by equations (|) and (|TI]) in the 



finite-temperature Hartree(-Fock) approximation, which is similar to the finite temperature 
local density approximation actually used in our numerical calculations. In a transparent, 
short notation, we describe the system by a second-quantized Hamiltonian 7i containing the 
kinetic energy operators tj of the particles present in the plasma, their two-body interactions 
Vij, and the additional one-body operators Wi = —Z/vi or +ZjZ/ri representing the 
contributions of the additional nuclear charge Z at the center. The temperature T enters 



the formalism via the the usual Boltzman factor exp{—f37i) where (3 = 1/T. Actually, 
because the average electronic density and background positive charge density are fixed 
parameters of the problem, Ti must be replaced by the usual constrained Hamiltonian. 
This will read, in a short notation, 7i — ^'M . Here /i' is a several- component chemical 
potential because of the several components in the plasma (i.e., electrons and various 
background ions). Accordingly the particle number operator M must be understood as a 
several-component operator. For simplicity, however, the following equations deal with the 
electrons only. This is because the positive ions in the background may consist of bosons 
as well as fermions, and we want to spare the reader the cumbersome symmetrization or 
antisymmetrization formulae for such backgrounds. In any case, the present paper describes 
the background density classically, and this density is not high enough to demand exchange 
terms. 

The independent particle ansatz approximates the eigenstates oi 7i — fi'Af by 
antisymmetrized products of single particle states (orbitals) \i), with an approximate 
spectrum made of sums of independent particle levels rji. The infinite set of orbitals \i) and 
eigenvalues rji define a one-body operator Hq — yu'A/", presumed optimal from the point of 
view of a variational principle: i.e., stationarity of the grand potential or free energy F. 

In second quantization, the many-body density matrix used as a trial density operator 
(noted in the following by quantities carrying a subscript t) is 

^ exp[-f3{no- fiW)] ^ exp{-(3J2iViclci) 

where q and q are the usual fermionic creation and annihilihation operators for orbital \i). 
For electrons the normalization denominator Zt is 

Zi = n[l + exp(-/377,)], (16) 
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so the trial density becomes 

_ -n- exp(-pr]iclci) 

The true density matrix D = exp[—(3{7i — yu'A/')]/Tr exp[—f3{7i — /u'A/")], minimizes the true 
grand potential, fi = Tr DiTi — ^'M + j3~^ InZ^). However, with an independent particle 
approximation one must be content with minimizing 

F = Tr A(H-/i'A/' + /3-MnA). (18) 

It is easy to show that 

Ti Dtc\cic]cj = fifj, if i^i, (19) 

where the Fermi occupation numbers are fi = Tr Dtclci = 1/[1 + exp{Pr]i)]. Furthermore 

Tr Dt {-f37],cU - ln[l + exp(-/5r/,)]} = /, In /. + (1 - /.) ln(l - /,). (20) 

Finally, because H = {T + W) + V is the sum of a one-body operator T + W = J2iti + J2iWi 
and a two-body operator V = J2i>j Vij, we obtain 

i i,j i 

Here the tilde v means that the matrix element {ij\v\ij) is antisymmetrized. 
The functional derivative of F with respect to {i\ is then 

^ = /. (t + w-fi' + Y.f.i.Mm, (22) 

where one recognizes the action of a mean field potential, including both direct and 
exchange terms 

i(\t)^Y.fA-M^j)- (23) 

Stationarity of F with respect to \i) then gives the "finite temperature Hartree-Fock" 
equations, 

{t + w-fi' + U)\i)=r]i\i), Vi, (24) 
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where one recognizes that rji is the Lagrange multipher for the normahzation of the orbital. 

In the same way, the derivative of F with respect to rji, or as well /», yields 

dF 

gJ^=V^ + r' In [/./I - /.] , (25) 

which vanishes [see Eq. (20)]. It can thus be concluded that 

F = -/3-MnZi-TrAV, (26) 

(which differs from —(3~^\n.Zt) is stationary with respect to variations of Dt. Therefore, 

since W is proportional to the additional charge Z at the center, the derivative of F with 

respect to Z must be given by 

dF W 

Notice also that W is a local potential. Hence, in the coordinate representation, only the 
diagonal part of Dt (the classical density) is needed to calculate dF/dZ, 

^^ ■ ^ '"W(0-Pe(rOl, (28) 



47r / r dr' 
dZ 



which is the differential form of Eqs. (11-12). Here Pe{r') and Pp{r') are diagonal matrix 



elements {r'\Dt\r') in the electron and positive background sectors, respectively. It will 
be noticed that Eq. (28) gives the background density as well as the electronic density, 
while the preceding equations, Eqs. (15-27), accounted for the electrons only. In view of 
the simplicity of the result due to the electrons, this reinstatement of the background 
contribution is trivial. Notice also that for the one-center problem, Dt is rotationally 
invariant, hence no vector label f' is needed. 

The result, Eq. (28) is simply a reformulation of Eqs. (11-12). The method of Matsubara 
poles used by |Watrous et al. 1999| via the finite temperature one-body Green's function, is 
perfectly suited to this local calculation of the density and avoids an explicit solution of 
Eqs. (24). Indeed, as discussed in detail by [Watrous et al. 199^ , a local density such as 
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Pe{f') can be directly derived from diagonal matrix elements (r'|(?7 — t~w — fi' — W)^^|r') of 
the one-body Green's function. Such matrix elements are integrated on a suitable contour 
in the ?7-complex plane. Once dF/dZ is known for all values of Z smaller than Zi + Z2, 
trivial integrals provide the screening f/o according to Eq. (10). 



5. NUMERICAL RESULTS FOR THE SOLAR CORE 

Shielding calculations have been performed for Z = 1,2,4,6,8 at a density and 
temperature relevant to energy and neutrino production in the sun. From these calculations 
the fusion enhancement factors are calculated and compared with the results of other 
researchers. IRicci et al. 19951 have shown that enhancement factors for the relevant fusion 



reactions are insensitive to the location within the solar core. We use here the set of 
parameters similar to those employed by [Bahcall and Pinsonneault 1995| , corresponding to 
R/Rq = 0.06. In atomic units, they are T = 47, n^ = 7.63, X = 0.432, Y = 0.568. 
The results are displayed in Table 1 for (p{Z, r = 0) and Ape{Z, r = 0) = Pe{Z, r = 0) — Ue- 
These quantities for other Z-values can be obtained by interpolation. (j){Z, r = 0) was fitted 
to a fourth order polynomial. 



4 



0(Z,r = O)^-^c„Z". (29) 

n=l 

The integral of this quantity then yields the free energy, 

F(Z)^-^c„Z"+V^, (30) 



n=l 



with ci = 2.0275, Ca = -0.03661, C3 = 0.00594, C4 = -0.000103. Note that (j){Z,r = 0) is 
nearly linear in Z. 

Table II gives the fusion enhancement factors for several reactions of interest in the 
solar neutrino problem. The deviation from unity is due to shielding. Note that the various 
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calculations exhibit a scatter of 10 to 30% and more in the deviation from unity. Although 
our factors are rather close to the recent calculations of [Gruzinov and Bahcall 19"98| , they 
are usually somewhat lower. 



6. Summary and Conclusions 



The finite temperature Green's function method of |Watrous et al. 1999| has been 



applied to the problem of screening of the nucleus-nucleus interaction in the solar plasma. 
The method is based on a self consistent, finite temperature Kohn-Sham formalism with 
the local density approximation for exchange and correlation. Atomic bound states are 
included on an equal footing with continuum states. Fusion enhancement factors (over 
pure Coulomb) are calculated for various relevant nuclear reactions at mean conditions 
(temperature and density) of the solar core. Comparisons with several other calculations 
are presented. 

The method appears to have no restriction with respect to temperature or density for 
stars in "non catastrophic," thermal equilibrium states. When implementing calculations 
for other systems, one may need to extend the range of parameters for the local 
exchange-correlation term, which is quite small for our solar system. 
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Table 1: Results of solar screening calculations for various values of Z. 
calculated values; the free energy F is based on the polynomial fit. 



1 


2.007 


2.791 


-1.008 


2 


3.944 


6.829 


-4.009 


3 


— 


— 


-8.962 


4 


7.884 


20.76 


-15.82 


5 


— 


— 


-24.54 


6 


12.00 


49.77 


-35.14 


7 


— 


— 


-47.65 


8 


16.51 


108.6 


-62.19 



and Ape are 



Table 2: Fusion enhancement factors for the solar interior (This Work), compared with other 
screening calculations: WES is Salpeter's weak screening approximation (|Salpeter 1954 ); Mit 
is Mitler's Thomas-Fermi-like modelf ^vTitTer 1977|) ; GDGC is due to [Graboske et al. 1973 



CSK is the dynamic screening model of [Uarraro et al. 1988| ; GB is due to [Gruzinov and 
Bahcall 1998"| . 



reaction This Work WES Mit GDGC CSK GB 

1.049 1.045 1.049 1.038 1.053 

1.213 1.176 1.115 1.158 1.224 

1.213 1.171 1.112 1.169 1.166 

1.402 1.293 1.191 1.324 1.393 



p + p 


1.043 


ie+He 


1.181 


Be+p 


1.183 


N+p 


1.356 
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